Load libraries

library(tidyverse) # data manipulation
library(ggpubr) # producing data exploratory plots
library(modelsummary) # descriptive data 
library(glmmTMB) # running generalised mixed models 
library(DHARMa) # model diagnostics 
library(performance) # model diagnostics  
library(ggeffects) # partial effect plots 
library(car) # running Anova on model 
library(emmeans) # post-hoc analysis

Import data

m1 <- read_csv("import_data/1_month_size_data_2022_2023.csv") |> 
  mutate(across(1:15,factor)) |> 
  mutate(STANDARD_LENGTH =LENGTH, 
         .keep = "unused") |> 
  select(!(NOTES)) |> 
  select(1:15,"STANDARD_LENGTH","MASS")
           
m2 <- read_csv("import_data/2_month_size_data_2022_2023.csv") |> 
  mutate(across(1:15,factor)) |> 
  mutate(STANDARD_LENGTH =LENGTH, 
         .keep = "unused") |> 
  select(!(NOTES)) |> 
  select(1:15,"STANDARD_LENGTH","MASS")
m2.5 <- read_csv("import_data/2-5_month_size_data_2022_2023.csv") |> 
  mutate(across(1:15,factor)) |> 
  mutate(STANDARD_LENGTH =LENGTH, 
         .keep = "unused") |> 
  select(!(NOTES)) |> 
  select(1:15,"STANDARD_LENGTH","MASS")
adult <- read_csv("import_data/adult_size_2022_2023.csv") |> 
  mutate(across(1:3,factor), 
         MALE = FISH_ID, 
         FEMALE = FISH_ID, 
         POPULATION = str_sub(FISH_ID, 2,4), 
         POPULATION = case_when(POPULATION == "ARL" ~ "Arlington Reef", 
                                POPULATION == "SUD" ~ "Sudbury Reef",
                                POPULATION == "VLA" ~ "Vlassof cay",
                                POPULATION == "PRE" ~ "Pretty patches", 
                                TRUE ~ POPULATION)) |> 
  left_join(select(m1, c("MALE","TEMPERATURE")), 
             by="MALE") |> 
  left_join(select(m1, c("FEMALE","TEMPERATURE")), 
             by="FEMALE") |>
  distinct() |> 
  mutate(TEMPERATURE = coalesce(TEMPERATURE.x, TEMPERATURE.y)) |> 
  drop_na(TEMPERATURE) |> 
  select(-c("TEMPERATURE.x","TEMPERATURE.y"))

Data manipulation

m1_df <- m1 |> 
  left_join(select(adult, c("MALE", "SL", "MASS")), 
            by ="MALE") |> 
  mutate(SL_MALE =SL, 
         MASS_MALE =MASS.y, 
         .keep = "unused") |>
  left_join(select(adult, c("FEMALE", "SL", "MASS")), 
            by ="FEMALE") |> 
  mutate(SL_FEMALE =SL, 
         MASS_FEMALE =MASS, 
         .keep ="unused") |> 
  mutate(SL_MIDPOINT = (SL_MALE+SL_FEMALE)/2, 
         MASS_MIDPOINT = (MASS_MALE+MASS_FEMALE)/2)

m2_df <- m2 |> 
  left_join(select(adult, c("MALE", "SL", "MASS")), 
            by ="MALE") |> 
  mutate(SL_MALE =SL, 
         MASS_MALE =MASS.y, 
         .keep = "unused") |>
  left_join(select(adult, c("FEMALE", "SL", "MASS")), 
            by ="FEMALE") |> 
  mutate(SL_FEMALE =SL, 
         MASS_FEMALE =MASS, 
         .keep ="unused") |> 
  mutate(SL_MIDPOINT = (SL_MALE+SL_FEMALE)/2, 
         MASS_MIDPOINT = (MASS_MALE+MASS_FEMALE)/2)

m2.5_df <- m2.5 |> 
  left_join(select(adult, c("MALE", "SL", "MASS")), 
            by ="MALE") |> 
  mutate(SL_MALE =SL, 
         MASS_MALE =MASS.y, 
         .keep = "unused") |>
  left_join(select(adult, c("FEMALE", "SL", "MASS")), 
            by ="FEMALE") |> 
  mutate(SL_FEMALE =SL, 
         MASS_FEMALE =MASS, 
         .keep ="unused") |> 
  mutate(SL_MIDPOINT = (SL_MALE+SL_FEMALE)/2, 
         MASS_MIDPOINT = (MASS_MALE+MASS_FEMALE)/2)

Exploratory data analysis

STANDARD_LENGTH

plot1 <- ggplot(m1_df, aes(x=SL_MALE, y=STANDARD_LENGTH, color=TEMPERATURE)) +
  geom_point(alpha=0.05) + 
  stat_smooth(method = "lm") + 
  theme_classic()

plot2 <- ggplot(m1_df, aes(x=SL_FEMALE, y=STANDARD_LENGTH, color=TEMPERATURE)) +
  geom_point(alpha=0.05) + 
  stat_smooth(method = "lm") + 
  theme_classic()

plot3 <- ggplot(m1_df, aes(x=SL_MIDPOINT, y=STANDARD_LENGTH, color=TEMPERATURE)) +
  geom_point(alpha=0.05) + 
  stat_smooth(method = "lm") + 
  theme_classic()

ggarrange(plot1, plot2, plot3, 
          nrow =1, 
          ncol =3, 
          common.legend = TRUE)

MASS

plot1 <- ggplot(m1_df, aes(x=MASS_MALE, y=MASS.x, color=TEMPERATURE)) +
  geom_point(alpha=0.05) + 
  stat_smooth(method = "lm") +
  ylim(0,0.15) +
  theme_classic()

plot2 <- ggplot(m1_df, aes(x=MASS_FEMALE, y=MASS.x, color=TEMPERATURE)) +
  geom_point(alpha=0.05) + 
  stat_smooth(method = "lm") + 
  ylim(0,0.15) +
  theme_classic()

plot3 <- ggplot(m1_df, aes(x=MASS_MIDPOINT, y=MASS.x, color=TEMPERATURE)) +
  geom_point(alpha=0.05) + 
  stat_smooth(method = "lm") + 
  ylim(0,0.15) +
  theme_classic()

ggarrange(plot1, plot2, plot3, 
          nrow =1, 
          ncol =3, 
          common.legend = TRUE)

Descriptive statistics

counts

Adults

datasummary(Factor(POPULATION) ~ Factor(TEMPERATURE), 
            data=adult, 
            fmt = "%.0f")
tinytable_fp19oyrgwc7ht7i43s5t
POPULATION 27 28.5 30
Arlington Reef 8 8 8
Pretty patches 4 6 6
Sudbury Reef 4 4 2
Vlassof cay 6 2 6

1-months

datasummary(Factor(POPULATION) ~ Factor(TEMPERATURE), 
            data=m1_df, 
            fmt = "%.0f")
tinytable_jp9kynvxbq8irggnsd6s
POPULATION 27 28.5 30
Arlington Reef 231 105 202
Pretty Patches 116 82 142
Sudbury Reef 117 55 47
Vlassof Cay 114 60 120

2-months

datasummary(Factor(POPULATION) ~ Factor(TEMPERATURE), 
            data=m2_df, 
            fmt = "%.0f")
tinytable_170h6aqgtgu96rziuc2g
POPULATION 27 28.5 30
Arlington Reef 198 108 152
Pretty Patches 113 83 134
Sudbury Reef 113 60 56
Vlassof Cay 77 58 113

2.5-months

datasummary(Factor(POPULATION) ~ Factor(TEMPERATURE), 
            data=m2.5_df, 
            fmt = "%.0f")
tinytable_8rzwk08isu9rif8eh67h
POPULATION 27 28.5 30
Arlington Reef 239 124 211
Pretty Patches 100 83 133
Sudbury Reef 111 53 50
Vlassof Cay 102 60 106

standard length

Adults

datasummary(Factor(TEMPERATURE) ~ SL * (NUnique + mean + median + min + max + sd + Histogram), 
            data = drop_na(adult, SL),  
            fmt = "%.2f")
tinytable_qisr8skk406du9q3ptgv
TEMPERATURE NUnique mean median min max sd Histogram
27 21 97.14 100.60 84.76 104.42 7.00 ▃▂▁▁▂▁▆▇
28.5 20 91.72 93.83 72.47 100.00 7.57 ▁▁▂▃▅▇▃
30 22 92.52 93.22 80.34 101.22 6.44 ▂▅▂▃▂▃▇▂▇▅

1-months

datasummary(Factor(TEMPERATURE) ~ STANDARD_LENGTH * (NUnique + mean + median + min + max + sd + Histogram), 
            data = drop_na(m1_df, STANDARD_LENGTH),  
            fmt = "%.2f")
tinytable_9jhrlvmjoqnmpogsgu9i
TEMPERATURE NUnique mean median min max sd Histogram
27 402 12.94 12.99 8.25 18.05 1.86 ▁▂▃▆▇▇▄▃▁
28.5 220 13.26 13.23 8.74 17.87 1.45 ▂▄▇▅▆▁
30 344 13.51 13.48 7.73 17.76 1.63 ▁▃▅▇▆▄▂

2-months

datasummary(Factor(TEMPERATURE) ~ STANDARD_LENGTH * (NUnique + mean + median + min + max + sd + Histogram), 
            data = drop_na(m2_df, STANDARD_LENGTH),  
            fmt = "%.2f")
tinytable_zd2zc1ps9q25huiiyldf
TEMPERATURE NUnique mean median min max sd Histogram
27 368 20.36 20.81 10.73 27.53 2.41 ▁▁▂▆▇▃
28.5 264 20.27 20.47 12.87 26.83 2.40 ▁▁▄▅▇▆▂▁
30 351 20.08 20.28 11.69 26.14 2.58 ▁▁▁▁▄▇▇▄▂

2.5-months

datasummary(Factor(TEMPERATURE) ~ STANDARD_LENGTH * (NUnique + mean + median + min + max + sd + Histogram), 
            data = drop_na(m2.5_df, STANDARD_LENGTH),  
            fmt = "%.2f")
tinytable_pm1dtouaw4oohq90yey1
TEMPERATURE NUnique mean median min max sd Histogram
27 439 24.15 24.56 11.92 31.70 3.02 ▁▁▂▄▇▅▂
28.5 283 23.65 24.05 14.69 30.15 3.00 ▁▁▁▃▅▅▇▆▃▁
30 382 23.58 23.94 13.74 31.04 2.61 ▁▁▂▄▇▇▃▁

mass

Adults

datasummary(Factor(TEMPERATURE) ~ MASS * (NUnique + mean + median + min + max + sd + Histogram), 
            data = drop_na(adult, MASS),  
            fmt = "%.2f")
tinytable_qr4c9h06i2qtkm376kr7
TEMPERATURE NUnique mean median min max sd Histogram
27 21 46.93 50.63 29.85 59.93 10.30 ▇▁▃▄▄▆▄
28.5 20 38.81 42.36 16.28 53.26 10.26 ▁▁▄▃▄▇▆▁
30 22 39.94 39.62 23.91 57.31 9.43 ▅▇▂▇▇▇▂▇▅▂

1-months

datasummary(Factor(TEMPERATURE) ~ MASS.x * (NUnique + mean + median + min + max + sd + Histogram), 
            data = drop_na(m1_df, MASS.x),  
            fmt = "%.2f")
tinytable_7tdoe1s0wasz4pfaingw
TEMPERATURE NUnique mean median min max sd Histogram
27 437 0.06 0.06 0.01 0.22 0.03 ▃▇▇▄▂▁
28.5 248 0.07 0.06 0.02 0.19 0.02 ▁▅▇▄▃
30 389 0.07 0.07 0.01 0.15 0.02 ▄▅▇▆▅▃▂▁

2-months

datasummary(Factor(TEMPERATURE) ~ MASS.x * (NUnique + mean + median + min + max + sd + Histogram), 
            data = drop_na(m2_df, MASS.x),  
            fmt = "%.2f")
tinytable_0wcs23i8tz9hnshau0vr
TEMPERATURE NUnique mean median min max sd Histogram
27 460 0.27 0.27 0.03 0.68 0.10 ▁▂▄▇▆▃▁
28.5 297 0.28 0.27 0.06 0.64 0.10 ▁▃▇▇▇▄▁▁▁
30 425 0.28 0.27 0.04 0.62 0.11 ▂▂▅▇▇▄▂▁▁

2.5-months

datasummary(Factor(TEMPERATURE) ~ MASS.x * (NUnique + mean + median + min + max + sd + Histogram), 
            data = drop_na(m2.5_df, MASS.x),  
            fmt = "%.2f")
tinytable_ysf3my51ycjdgj0kh11m
TEMPERATURE NUnique mean median min max sd Histogram
27 529 0.46 0.46 0.04 1.00 0.17 ▁▂▄▇▇▆▄▁▁
28.5 314 0.45 0.44 0.08 1.02 0.18 ▂▄▇▇▆▅▃▁▁
30 477 0.45 0.45 0.08 0.96 0.15 ▁▂▄▆▇▄▂▁